## Attainment Status Analysis (Results are estimated by lfe package version 2.8-7)
## These results are reported in Table 3
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(readxl)
library('usmap')
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
## Generate county info
data3$StateName=as.character(data3$StateName)
data3$County=as.character(data3$County)
county_list=unique(data3[c("StateName","County")])
## Fix county names to match non-attainment data
county_list$County1=county_list$County
county_list=subset(county_list, !(is.na(County1)))
county_list$County1[which(county_list$County1=="Miami Dade")]="Miami-Dade"
county_list$County1[which(county_list$County1=="St Lucie")]="St. Lucie"
county_list$County1[which(county_list$StateName=="Louisiana")]=paste(county_list$County1[which(county_list$StateName=="Louisiana")], ' Parish', sep='')
county_list$County1[which(county_list$County1=="St Mary Parish")]="St. Mary Parish"
county_list$County1[which(county_list$County1=="DeSoto Parish")]="De soto Parish"
county_list$County1[which(county_list$County1=="St James Parish")]="St. James Parish"
county_list$County1[which(county_list$County1=="St Charles Parish")]="St. Charles Parish"
county_list$County1[which(county_list$County1=="St Bernard Parish")]="St. Bernard Parish"
county_list$County1[which(county_list$County1=="Queen Annes")]=paste("Queen Anne's")
county_list$County1[which(county_list$County1=="Prince Georges")]=paste("Prince George's")
county_list$County1[which(county_list$County1=="St Clair")]="St. Clair"
county_list$County1[which(county_list$County1=="St Francois")]="St. Francois"
county_list$County1[which(county_list$County1=="St Louis")]="St. Louis"
county_list$County1[which(county_list$County1=="St Charles")]="St. Charles"
county_list$County1[which(county_list$County1=="St Louis City")]="St. Louis City"
county_list$County1[which(county_list$County1=="DeWitt")]="De Witt"
county_list$County1[which(county_list$County1=="Dewitt")]="De Witt"
county_list$County1[which(county_list$County1=="St Joseph")]="St. Joseph"
county_list$County1[which(county_list$County1=="St Croix")]="St. Croix"
county_list$County1[which(county_list$County1=="St Lawrence")]="St. Lawrence"
county_list$fips=NA
for (i in 1:nrow(county_list)) {
  county_list$fips[i]=fips(county_list$StateName[i], county_list$County1[i])
  print(i)
}
county_list=county_list[c("StateName","County","fips")]

data3=left_join(data3, county_list, by=c("StateName","County"))
data3$fips=as.numeric(data3$fips)

########################
## Extract attainment information
nonattain=read_xls("raw data/nonattain.xls") ## From EPA NAAQS nonattainment status
nonattain$fips=as.numeric(nonattain$fips_state)*1000+as.numeric(nonattain$fips_cnty)
nonattain=nonattain[c('fips',"pollutant","pw_2013","pw_2014","pw_2015","pw_2016","pw_2017","pw_2018")]
nonattain$Year_2013="Attain"
nonattain$Year_2013[which(nonattain$pw_2013=='P')]='Partial'
nonattain$Year_2013[which(nonattain$pw_2013=='W')]='Whole'
nonattain$Year_2014="Attain"
nonattain$Year_2014[which(nonattain$pw_2014=='P')]='Partial'
nonattain$Year_2014[which(nonattain$pw_2014=='W')]='Whole'
nonattain$Year_2015="Attain"
nonattain$Year_2015[which(nonattain$pw_2015=='P')]='Partial'
nonattain$Year_2015[which(nonattain$pw_2015=='W')]='Whole'
nonattain$Year_2016="Attain"
nonattain$Year_2016[which(nonattain$pw_2016=='P')]='Partial'
nonattain$Year_2016[which(nonattain$pw_2016=='W')]='Whole'
nonattain$Year_2017="Attain"
nonattain$Year_2017[which(nonattain$pw_2017=='P')]='Partial'
nonattain$Year_2017[which(nonattain$pw_2017=='W')]='Whole'
nonattain$Year_2018="Attain"
nonattain$Year_2018[which(nonattain$pw_2018=='P')]='Partial'
nonattain$Year_2018[which(nonattain$pw_2018=='W')]='Whole'
## Generate nonattainment list
nonattain=subset(nonattain, Year_2018!="Attain"|Year_2017!="Attain"|Year_2016!="Attain"
                 |Year_2015!="Attain"|Year_2014!="Attain"|Year_2013!="Attain")
nonattain_all=subset(nonattain, Year_2018!="Attain"&Year_2017!="Attain"&Year_2016!="Attain"
                     &Year_2015!="Attain"&Year_2014!="Attain"&Year_2013!="Attain")
nonattain=unique(nonattain$fips)
nonattain_all=unique(nonattain_all$fips)
## Drop counties switching between attainment and nonattainment status
switcher=subset(nonattain, !(nonattain%in%nonattain_all))
data3$fips=as.numeric(data3$fips)
data3=subset(data3, !(fips%in%switcher))
# Define attainment dummy
data3$attain=1
data3$attain[which(data3$fips%in%nonattain)]=0


data3$shutdown_attain=0
data3$shutdown_attain[which(data3$shutdown==1&data3$attain==1)]=1

data3$shutdown_nonattain=0
data3$shutdown_nonattain[which(data3$shutdown==1&data3$attain==0)]=1

########################
# Generate sample without post-shutdown period
data4=subset(data3, T<99)
########################
# Generate the short sample
data5=subset(data4, T>40)

##############################################################################################################################
########################
## DID regressions
data_base1=data4
data_base2=data5
all_AOD_long=felm(AOD~ shutdown_attain + shutdown_nonattain  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data_base1)
summary(all_AOD_long)
all_SO2_long=felm(SO2_MASS..tons.~ shutdown_attain + shutdown_nonattain  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data_base1)
summary(all_SO2_long)
all_NOX_long=felm(NOX_MASS..tons.~ shutdown_attain + shutdown_nonattain  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data_base1)
summary(all_NOX_long)
all_AOD_short=felm(AOD~ shutdown_attain + shutdown_nonattain  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs.  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data_base2)
summary(all_AOD_short)
all_SO2_short=felm(SO2_MASS..tons.~ shutdown_attain + shutdown_nonattain + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data_base2)
summary(all_SO2_short)
all_NOX_short=felm(NOX_MASS..tons.~shutdown_attain + shutdown_nonattain  + ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                   data = data_base2)
summary(all_NOX_short)
## Output results
stargazer(all_SO2_long, all_NOX_long,all_AOD_long,all_SO2_short, all_NOX_short, all_AOD_short)






